home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_11_06 / 1106032a < prev    next >
Text File  |  1993-02-15  |  21KB  |  658 lines
  1. /******************************************************
  2.  * RPFT.C -- Curve fitting with optional extrapolation.
  3.  * Fits either polynomial or rational function. Based
  4.  * on algoritms defined in Press, et al., "Numerical
  5.  * Recipes", 2nd ed., Cambridge University Press, 1992
  6.  * Developed using the Borland C compiler.
  7.  *
  8.  *   Lowell Smith
  9.  *   3368 Crestwood Drive
  10.  *   Salt Lake City, UT  84109-3202
  11.  *****************************************************/
  12. #include <stdio.h>
  13. #include <math.h>
  14. #include <string.h>
  15. #include <dir.h>
  16. #include <ctype.h>
  17. #include <stddef.h>
  18. #include <stdlib.h>
  19.  
  20. #define NPFAC 8
  21. #define PI 3.141592653589793
  22. #define PIO2 (PI/2.0)
  23. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  24. void ratlsq(double xs[],double fs[],int npt,int mm,
  25.       int kk, double cof[], double *dev);
  26. double *dvector(long nl, long nh, const char *ner);
  27. double **dmatrix(long nrl,long nrh,long ncl,long nch,
  28.          const char *ner);
  29. void free_dvector(double *v, long ncl);
  30. void free_dmatrix(double **m, long nrl, long ncl);
  31.  
  32. /*****************************************************/
  33. void main(int argc, char **argv)
  34. { double a,b,*xs,*fs,cof[42],dev=0.,sumd,x,yy;
  35.   int i,j,nn,nd,ncof,dig,npt=0,xln,yln;
  36.   char nar[90];
  37.   char filenam[80];
  38.   void commd(int argc, char**argv, double *a,double *b,
  39.         int *nn, int *nd, int *dig, int *xln,
  40.         int *yln, char *filenam);
  41.   void inpt(int nn,int npt,double a, double b, int xln,
  42.         int yln, double *xs,double *fs,char filenam[]);
  43.   void pcshft(double a,double b,double d[],int n);
  44.   void chebpc(double c[],double d[],int n);
  45.  
  46.   commd(argc,argv,&a,&b,&nn,&nd,&dig,&xln,&yln,
  47.         filenam);
  48.   if (xln) { a=log(a);  b=log(b); }
  49.   if (nd)            /* Fit rational function */
  50.   { ncof=nn+nd+1;   npt=NPFAC*ncof;
  51.     xs=dvector(1L,(long)npt,"xs in main");
  52.     fs=dvector(1L,(long)npt,"fs in main");
  53.     inpt(0, npt, a, b, xln, yln, xs, fs, filenam);
  54.     ratlsq(xs,fs,npt,nn,nd,&cof[0],&dev);
  55.     fprintf(stderr,"Est. max error = %.7G\n",dev);
  56.     free_dvector(xs,1L);
  57.     free_dvector(fs,1L);
  58.   }
  59.   else               /* Fit polynomial */
  60.   { ncof=nn+5; xs=dvector(0L,(long)ncof,"xs in main");
  61.     fs=dvector(0L,(long)ncof,"fs in main");
  62.     inpt(ncof, npt, a, b, xln, yln, xs, fs, filenam);
  63.     chebpc(fs,cof,nn+1);
  64.     pcshft(a,b,cof,nn+1);
  65.     free_dvector(xs,0L); free_dvector(fs,0L);
  66.   }
  67.   for (i=0;i<=nn+nd;++i)
  68.   { sprintf(nar,"%.*G ",dig,cof[i]);
  69.     sscanf(nar,"%lf",&cof[i]);
  70.   }
  71.   if (nd) printf("\n      Rational function coefficien"
  72.      "ts\n   Numerator               Denominator\n\n");
  73.   else printf("\n  Polynomial\n  Coefficients\n\n");
  74.   for (i=0;i<max(nd,nn)+1;++i)
  75.   { if (i<=nn) printf("%-24.16G",cof[i]);
  76.     else printf("%*s",24," ");
  77.     if (i<nd) printf("%-24.16G\n",cof[i+nn+1]);
  78.     else printf("\n");
  79.   }
  80.   for(;;)            /* Calculate trial values */
  81.   { i=-1;
  82.     fprintf(stderr,
  83.           "Enter E to quit or a trial X value: ");
  84.     i=scanf("%lf",&x);
  85.     if (i<1) exit(1); if (xln) x=log(x);
  86.     yy=cof[nn]; for (j=nn-1;j>=0;j--) yy=yy*x+cof[j];
  87.     if (nd)
  88.     { for (sumd=0.,j=nn+nd;j>=nn+1;j--)
  89.                  sumd=(sumd+cof[j])*x;
  90.       yy=yy/(1.0+sumd);
  91.     }
  92.     if (yln) yy=exp(yy); if (xln) x=exp(x);
  93.     fprintf(stderr,"For X=%.8G  Y=%.8G\n",x,yy);
  94.   }
  95. }
  96. /******************************************************
  97.  * COMMD - Parses the command line
  98.  *****************************************************/
  99. void commd(int argc, char**argv, double *a,double *b,
  100.         int *nn, int *nd, int *dig, int *xln,
  101.         int *yln, char *filenam)
  102. { struct ffblk ffblk;
  103.   int i,j,k,l;
  104.   void help(char *msg);
  105.  
  106.   *a=1.11e300; *b=-1.11e300; *nn=-32768; *nd=-32768;
  107.   *dig=7, *xln=0, *yln=0;
  108.   if (argc<2) help(""); filenam[0]=0;
  109.   for (i=1,k=0;i<argc;k=0,++i)
  110.   { for (j=0; j<(l=strlen(argv[i])); ++j)
  111.     { argv[i][j]=toupper(argv[i][j]);
  112.       if (argv[i][j]=='=') ++k;
  113.     }
  114.     if (k)
  115.     { if (k>1 || (!isdigit(argv[i][l-1])
  116.                       && argv[i][l-1] !='.'))
  117.       { fprintf(stderr,"Invalid command line parameter"
  118.                  " %s\n",argv[i]);
  119.         help("");
  120.       }
  121.       if (!strncmp(argv[i],"LL=",3))
  122.       { k=sscanf(&argv[i][3],"%lf",a);
  123.         if (k<1) help("Invalid value for command line"
  124.                    " parameter LL\n");  continue;
  125.       }
  126.       if (!strncmp(argv[i],"UL=",3))
  127.       { k=sscanf(&argv[i][3],"%lf",b);
  128.         if (k<1) help("Invalid value for command line"
  129.                    " parameter UL\n");  continue;
  130.       }
  131.       if (!strncmp(argv[i],"ND=",3))
  132.       { k=sscanf(&argv[i][3],"%d",nd);
  133.         if (k<1) help("Invalid value for command line"
  134.                    " parameter ND\n");  continue;
  135.       }
  136.       if (!strncmp(argv[i],"NN=",3))
  137.       { k=sscanf(&argv[i][3],"%d",nn);
  138.         if (k<1) help("Invalid value for command line"
  139.                    " parameter NN\n");  continue;
  140.       }
  141.       if (!strncmp(&argv[i][0],"-DIG=",5))
  142.       { k=sscanf(&argv[i][6],"%d",dig);
  143.         if (k<1) help("Invalid value for optional"
  144.                    " command line parameter DIG\n");
  145.         continue;
  146.       }
  147.       fprintf(stderr,"Invalid command line parameter "
  148.                    "%s\n",argv[i]);
  149.       help("");
  150.     }
  151.     else
  152.     { if (!strncmp(&argv[i][0],"-XLN",4))
  153.       { *xln=1; continue; }
  154.       if (!strncmp(&argv[i][0],"-YLN",4))
  155.       { *yln=1; continue; }
  156.       if (!filenam[0] && argv[i][0] != '-')
  157.       { strcpy(filenam,argv[i]);
  158.         if (findfirst(filenam,&ffblk,0))
  159.         { fprintf(stderr,"Input file %s doesn't "
  160.             "exist\n",filenam); help("");
  161.         }
  162.       }
  163.       else
  164.       { fprintf(stderr,"Invalid command line parameter"
  165.                    " %s\n",argv[i]);   help("");
  166.       }
  167.     }
  168.   }
  169.   k=0;
  170.   if (*a>1.1e300)
  171.   { fprintf(stderr,"Parameter LL undefined\n"); ++k; }
  172.   if (*b<-1.1e300)
  173.   { fprintf(stderr,"Parameter UL undefined\n"); ++k; }
  174.   if (*nn==-32768)
  175.   { fprintf(stderr,"Parameter NN undefined\n"); ++k;
  176.     *nn=5;
  177.   }
  178.   if (*nd==-32768)
  179.   { fprintf(stderr,"Parameter ND undefined\n"); ++k;
  180.     *nd=5;
  181.   }
  182.   if (filenam[0]==0)
  183.     { fprintf(stderr,"Input file is not defined\n");
  184.       ++k;
  185.     }
  186.   if (*nn<1 || (*nd==0 && *nn>20) || (*nd && *nn>10))
  187.   { fprintf(stderr,"Invalid value %d for "
  188.             "command line parameter NN\n",*nn); ++k;
  189.   }
  190.   if (*nd<0 || *nd>10)
  191.   { fprintf(stderr,"Invalid value %d for "
  192.            "command line parameter ND\n",*nd); ++k;
  193.   }
  194.   if (*a>*b)
  195.   { fprintf(stderr,"Lower limit > upper limit\n");
  196.     ++k;
  197.   }
  198.   if (*dig<1 || *dig>20)
  199.   { fprintf(stderr,"Invalid value %d for "
  200.           "command line option DIG\n",*dig); ++k;
  201.   }
  202.   if (*xln && *a<=0.)
  203.   { fprintf(stderr,"Lower limit is <=0  with "
  204.         "logarithmic transform of X values\n"); ++k;
  205.   }
  206.   if (k) help("");
  207. }
  208. /******************************************************
  209.  * HELP - Prints the help screen
  210.  *****************************************************/
  211. void help(char *msg)
  212. { char *hlp[] =
  213.   { "USAGE: rpft LL=a UL=b NN=n ND=m [-DIG=n -XLN -YLN"
  214.     "] file","","   LL=a   a is the lower limit of the"
  215.     "region for fit","   UL=b   b is the upper limit o"
  216.     "f the region for fit","","   ND=m | m>0: rational"
  217.     " function, denominator degree m,","   NN=n |     "
  218.     " numerator degree n.  1<=m<=10  1<=n<=10","      "
  219.     "  | m=0: Polynomial degree n.  1<=n<=20",""," -DI"
  220.     "G=n   (optional) n = number of significant digits"
  221.     ,"              for coefficients on output. Defaul"
  222.     "t=7.","   -XLN   (optional) Transform X values to"
  223.     " ln(X)","   -YLN   (optional) Transform Y values "
  224.     "to ln(Y)","","   file   File with input X Y data "
  225.     "points, 1 per line","","All command line paramete"
  226.     "rs except DIG, XLN and YLN are ","required. They "
  227.     "may be in any order, upper or lower case."
  228.   };
  229.   int i,n=sizeof(hlp)/sizeof(char *)-1;
  230.  
  231.   fprintf(stderr,"%s\n",msg);
  232.   for (i=0; i<n; ++i) fprintf(stderr,"%s\n",hlp[i]);
  233.   fprintf(stderr,"%s",hlp[n]);  exit(1);
  234. }
  235. /******************************************************